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Abstract 

The Hamiltonian Mean Field (HMF) model is a prototype for systems with long- 
range interactions. It describes the motion of N particles moving on a ring, coupled 
through an infinite-range potential. The model has a second order phase transition 
at the energy density C/c = 3/4 and its dynamics is exactly described by the Vlasov 
equation in the N ^ oo limit. Its chaotic properties have been investigated in the 
past, but the determination of the scaling with N of the Lyapunov Spectrum (LS) 
of the model remains a challenging open problem. We here show that the N^^^'^ 
scaling of the Maximal Lyapunov Exponent (MLE), found in previous numerical and 
analytical studies, extends to the full LS; not only, scaling is "precocious" for the LS, 
meaning that it becomes manifest for a much smaller number of particles than the 
one needed to check the scaling for the MLE. Besides that, the A^~^/^ scaling appears 
to be valid not only for U > Uc, as suggested by theoretical approaches based on a 
random matrix approximation, but also below a threshold energy Ut ~ 0.2. Using 
a recently proposed method (GALI) devised to rapidly check the chaotic or regular 
nature of an orbit, we find that Ut is also the energy at which a sharp transition 
from weak to strong chaos is present in the phase-space of the model. Around this 
energy the phase of the vector order parameter of the model becomes strongly time 
dependent, inducing a significant untrapping of particles from a nonlinear resonance. 

Key words: Hamiltonian Mean Field (HMF) model, Lyapunov spectra, GALI 
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1 Introduction 



The dynamical and statistical properties of system s with long-range in teractions 



have recently attracted considerable attention (see ICampa et al.l (|2009r) for a re- 
cent review). The Hamiltonian Mean Field (HMF) model ( lAntoni fc Ruffol . Il995l ) is 
considered as a prototype for some of the observed dynamical effects. The model de- 
scribes the motion of fully coupled particles lying on a ring. The system presents 
a second order phase transition at the critical energy density Uc = Ec/N = 3/4 
from a low-energy phase, where the particles are clustered, to a high-energy gaseous 
phase where the particles are uniformly distributed on the ring. Several different 
effects have been studied for this model, which have been then found to be generic 
for a large class of systems with long-range interactions, including self-gravitating 
systems, unscreened plasmas, two-dimensional hydrodynamics, etc.. 



The question on which we concentrate in this paper is the behavior of the Lyapunov 
exponents. This question has been investigated for a long time, already in some of 
the first papers on the HMF model. However, a general understanding of the chaotic 
properties of this model remains a challenge. Since we will be mostly concerned with 
the — !■ oo behavior of the Lyapunov exponents, this question has a relevance also 
for the Vlasov equation, which is known to describe exactly the dynamics of the 
HMF model in this hmit flCampa et all . l2009h . 



Let us briefly summarize the state of the art knowledge on the Lyapunov expo- 
nents for the HMF model. The Maximal Lyapunov Exponent (MLE) increases with 
energy up to the critical energy density Uc and then decreases for larger energies 
( Yamaguchi . 19961). In the whole h igh energy phase it drops to zero in the large N 
limit as N~^^^ ( Latora et al.l . Il998l ). An analytical estimate of the MLE, whi ch takes 
into account microcano nical averages of suitable geometrical observables (jPettinil . 
2007 ). was proposed by Firpol ( 1998 ). A preliminary study of the scaling proper- 
ties of the Ly apunov Spectra (LS ) and of the Kolmogorov-Sinai (KS) entropy was 
performed by iLatora et al.l ( 1l999l ). Lyapun ov instability of a modified HMF model, 
the so-calle d g-HMF, was first s tudie d by lAnteneodo fc TsallisI ( 1l998l ) and further 
analyzed by lCampa et al.l ( 200lll2002l ). The A^~^/^ scaling has been shown to derive 



20011 : IValleios fc Anteneodol. 120021: lAnteneodo et al. 



from a random matrix approximat i on (iFirpo fc Ruffo . l200ll : lAnteneodo fc Vallejos 



20031 ). A supersymmetric ap- 



proach ( iTanase-Nicola fc Kurchanl . |2003|) suggests that the MLE might vanish in 
the N ^ oo for all energies. A tendency of the model towards integrability when N 
increases, meaning a v anishing MLE, has been recently emphasized using a Poincare 



return map approach (IBachelard et al.l . 120081 ). Lyapunov exponents and the corre- 



sponding eigenmodes of the HMF model have been also recently studied in the 
Vlasov A^ — oo limit ( jPaskauskas fc De Ninnol . |2009|) . 



The aim of this paper is to present a careful assessment of the scaling properties 
with A^ of the MLE, the LS and the KS entropy of the HMF model in various 
energy ranges and in conditions as close as possible to Boltzmann-Gibbs thermal 
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equilibrium. With respect to previous studies, we have put more effort in estabhshing 
the equihbrium state, kno wing that thi s is no t easily reached due to the presence 



of quasi-stationary states (jCampa et al.l . |2009[ ). These latter are out-of-equilibrium 



states in which the system remains trapped for a time that increases with A^, slowing 
down indefinitely the relaxation to thermal equilibrium. 

An alternative ap proach, designed for fast detection of chaotic orbits, has been 



recently proposed (jSkokos et al.l . 120071 ). The method, which has been called Gen- 
eralized Alignment Index (GALI), allows one to distinguish between regular and 
chaotic motion using short-time numerical integrations. The GALI method takes 
into account more than one deviation vector and measures the time evolution of 
the volume of the parallelepiped whose edges are these normalized vectors. In this 
paper we use the GALI method to detect the fraction of chaotic orbits on a constant 
energy surface of the HMF model. 

The paper is organized as follows. In Section [2] we briefly introduce the HMF model 
and we discuss the main features of the Boltzmann-Gibbs statistical equilibrium 
state. In Section [3] we recall how the LS is computed, in order to make the paper 
self-consistent. Afterwards, in Section |U we present the core results of this paper, by 
discussing the scaling with A^ of the MLE, the LS and the KS entropy. In Section[5]we 
describe the GALI chaos detection method and we test its efficiency by considering 
the HMF model with a small number of particles. In Section we use the GALI 
method for the estimation of the fraction of chaotic orbits in the model's phase 
space. Our conclusions are finally summarized in Section [71 



2 The HMF model and its statistical equilibrium 



The Hamiltonian Mean Field (HMF) model (lAntoni fc Ruffol Il995[ ) describes a sys- 
tem of point masses moving on a ring and interacting through an infinite-range 
potential. The Hamiltonian of the model is 



N 



N 



^ o 2N 



cos 



(1) 



1,3= 



where 6i G (tt, it] is the coordinate of the i-th particle and Pi its conjugate momen- 
tum. The system has only two constants of the motion: total energy H and total 
momentum J^iPi (iii the following we will always choose a zero total momentum 
without loosing generality). As the force vanishes when i = j, the particles do not 
collide but smoothly cross each other. The order parameter of this model is the 
"magnetization" vector 



1 



N 



^ = jj E(cos 9^, sin 9i) = (M„ My) , 

i=l 



(2) 
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and, after defining its phase = arctanMj^/Mc and modulus M = ^ + My, the 
equations of motion can be written as 



Oi = Pi 

(3) 

p. = -M^ sin Oi + My cos Oi = -Msin(6'i - 0) . 



We solve numerically the equations of motion of the HMF mod el using an optimized 
fourth-order symplectic integrator ( McLachlan fc Atela . 19921 ) with time step h = 



0.05, which typically gives a relative energy fiuctuation AE/E ~ 10" 
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The equ ilibrium statis t ical m echanics of the HMF model has been thoroughly stud- 



ied (see ICampa et al.l (120091 ) for a review). The model has a second order phase 
transition at the critical energy density Uc = 3/4 from a low-energy clustered phase 
where M 7^ to a high-energy gaseous phase where M = 0. Being this phase tran- 
sition of second order, it cannot be associated with a negati ye specific heat. Th e 



canonical and microcanonical ensemble give equivalent results (jCampa et al.l . l2009l ). 



In the mean- field limit, N 00, the dynamics of the HMF model is fully described 
by the single particle distribution function f{6,p,t), where {6,p) are the canoni- 
cally conjugate Eulerian coordinates of the single particle phase space. The single 
particle distribution function obeys a Vlasov equation. Among many stationary so- 
lutions of the Vlasov equation, a particular one is the Boltzmann-Gibbs equilibrium 
distribution 

where /3 = 1/T is the inverse temperature and Iq is the modified Bessel function of 
order zero. This distribution is also stable with M = for /3 < /3c = 2 and with 
M 7^ for /3 > Pc- Therefore, Pc can be interpreted as the inverse critical temperature 
of the second order phase transition, corresponding to the critical energy Uc- 

In order to reach the equilibrium distribution on a reasonably short time scale, 
we have initialized positions with a Gaussian distribution and we have computed the 
corresponding potential energy. We then determined the appropriate width of the 
Gaussian distribution of momenta which gave the energy we wanted to achieve. This 
initial state was then evolved for a long time in order for the equilibrium distribution 
of Eq. dl]) to be reproduced with a sufficient precision. For instance, for the = 100 
particle case, we can reproduce the equilibrium temperature with a relative error of 
a few percent, which is compatible with the statistical error, expected to be of the 
order N"^/"^. 
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3 The Lyapunov spectrum 



For the sake of completeness, we briefly rec all how the Lyapunov s pectrum of a 2N- 
dimensional flow is deflned and computed (IBenettin et all ll980al Jbl). with reference 
to the HMF model. 



The 2iV-dimensional phase-space coordinate is 

X = (6'i, . . .,9n,Pi, ■ ■ • ,P7v) 



(5) 



and the phase-space flow is generated by the system of autonomous flrst-order dif- 
ferential equations ([3]), which, using x, can be written as 



cix(t) 
dt 



F(x(t)) 



(6) 



where F is the velocity fleld of the flow. 

The evolution equations for the deviation vector 

w = (56*1, . . . , 69n, Spi, 6pn) 

are 



dw 



J(x(t))w 



where the 2N x 2N Jacobian matrix of the flow J(x(t)) = dF/d^ is given by 



J 




(7) 
(8) 

(9) 



with the N X N Jacobian Jij = —d'^V/dOidOj given by 



^ Jii = -cos{9,)M,-sm{9i)My + ^ ^^^^ 
^ Jij = cos(6'i - Oj), if i^j, 

and / the N x N identity matrix. 

The linear evolution equations are integrated with initial value w(0), which 
physically represents the "small" initial difference between two nearby orbits of 
Eq. (Q. Numerically, this integration is performed in parallel with that of the orbit 
x(t), since J depends on it. One typically uses the same integration algorithm. 

The Maximal Lyapunov Exponent (MLE) is then deflned as 

1 llwmii 

A(x(0)) = lim - In ) ( . (11) 
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Practically, since ||w(t)|| diverges exponentially, one has to renormalize it at reg- 
ular time intervals. The computation of the MLE is performed by averaging the 
renormalization factors. 

Furthermore, when varying ||w(0)||, one can obtain at most 2N Lyapunov exponents, 
which can be ordered in size 

Ai > A2... > A2^ , (12) 
and constitute the so-called Lyapunov Spectrum (LS), with Ai = A(x(0)). 
For Hamiltonian flows, the LS is endowed of the following symmetry property 

K = -\2N-i+l. 2 = l,2,...,2iV . (13) 

It is then enough to compute the first N exponents. Among them, for the HMF 
model, two are zero because of energy and momentum conservation. For a chaotic 
orbit all the others are typically positive. 

However, because of the numerical instability which tends to align the deviation 
vector w(t) along the most expanding direction, one always obtains Ai and has not 
access to the full LS. 



The tri ck to perform numerically the calculation of the full LS was found by lBenettin et al 
fll980al lbh and amounts to compute the hypervolume 



Vp(t) = wi(t) A W2(t) A ... A Wp(t), (14) 
for p deviation vectors wi(t), Wp(t). 
Then, the sum of the first p Lyapunov exponents 

A(p) = Ai + A2 + ... + Ap (15) 

is given by 



A(P)(x(0)) = lim-lnp^. (16) 
'^"^t l|Vp(0)|| ^ ' 

A recent review about the Lyapunov exponents and their calculation can be found 
in 



skokosi mm . 



Once the LS is known, the Kolmogorov-Sinai (KS) entropy, which is the rate of 
information production, is given by 

N 
i=l 
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4 Scaling with of the MLE, the LS and the KS entropy 



As mentioned in the Introduction, it has been proposed that in the whole high 
energy region U > Uc the MLE should vanish as A^~^/^. We want to show here that 
this scaling law is obeyed also by the full LS. 

In panels a) and c) of Fig. [T] we plot vs. i/N, increasing N from = 50 to 
A^ = 500, for U = 0.1 (panel a) and U = 1.2 (panel c). The total integration time 
was 0.5 X 10^ with time step h = 0.05, which provides a good convergence for the LS 
and a good accuracy for the energy conservation as well. The relative error in the 
determination of the MLE is of the order of 1%. We observe a general decreasing 
trend, showing a significant size dependence of the LS. This size dependence is 
almost completely eliminated when multiplying the LS by N^^^, as shown in panels 
b) and d) of the same figure. We observe a reasonably good "data collapse" in the 
upper and lower parts of the spectrum. However, the insets in panels b) and d) show 
that the size dependence is not yet completely eliminated and remains of the order 
of 10%. The insets reveal another interesting feature: while the convergence to an 
asymptotic LS is from above for the full LS at f/ = 0.1, at the energy U = 1.2 it 
is from above only for the largest exponents and it is instead from below for the 
smaller exponents. 






Fig. 1. LS vs. i/N for U = 0.1 (panel a) and U = 1.2 (panel c). For the same energies we 
also show the LS multiplied by N^^^ (panels b and d). The insets show details of the size 
dependence with more resolution. 

In Fig. [2] we show the A^ dependence of Ai for U = 0.1, which is well fitted by the 
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scaling law N~^^^, confirming the result obtained by looking at the full LS. However, 
one has to reach much larger system sizes (of the order of 10^) in order to check the 
scaling law, while the convergence to the asymptotic LS as increases is observed 
for much smaller system sizes (up to iV = 500). 



Othe r numerical experiments ( Latora et al. . 19981 : Anteneodo fc Tsallis . 1998 : Anteneodo et al. 



20031 ) show that this scaling law is present also for energies above Uc, but the quality 
of the data is never comparable with the one obtained here for U = 0.1. 



While the sca ling in the high energy range can be justified using a random matrix ap- 



proxi mation (ILatora et al.l . ll998l : lAnteneodo fc Tsallisl . ll998tlAnteneodo fc Valleios 



20011) and a theory using geometric properties of the phase-space (jPettinil . 12007 : 
Firpol Il998l ). no theoretical approach exists for the explanation of the N~^/^ scaling 
at low e nergie s . On t he other hand, the only theoretical approach valid in this energy 
region (iFirpd . Il998[ ) would predict a strictly positive MLE in the N ^ oo limit. 



U-0.1 


















■. 




v. 



Log(N) 

Fig. 2. Logarithm of the MLE vs. logA^ for U = 0.1, showing the A^~^/^ scaling. Both 
logarithms are in base 10. 

In Fig. [3^), we show the MLE as a function of U for various system sizes: = 
50, 100, 200, 500. One observes a sharp increase of the MLE around the energy Ut ~ 
0.2. We'll come back to comment about this energy value in Section |6l The MLE 
is peaked around the critical energy Uc and shows, in this range of sizes, a weak 
size dependence for U < Uc while for U > Uc it systematically and monotonically 
decreases with system size. In panel b) of the same figure, we check if the A^~^/^ 
scaling found for the LS works also for the MLE, by multiplying all MLE by the 
factor N^^^. It turns out that this scaling does not work well for this relatively small 
number of particles. An "indermediate" scaling with a factor N~^^^ looks slightly 
better in the high energy range, see Fig. [St. 

As a partial conclusion, we could state that the MLE shows a size dependence over 
the whole energy range, but the N~^^^ scaling does not emerge easily from the data. 
This is at variance with the analysis above of the LS, for which this scaling was 
more evident. This is why we speak of "precocious" scaling of the LS, meaning that 
the scaling is obtained for moderately large system sizes, while the scaling of the 
MLE, as shown e.g. in Fig. [21 needs much larger systems sizes to be numerically 



8 



emphasized. 




N-50 
N=100 
N=200 
N=500 



N-50 
N=100 
N=200 
N-50O 



Fig. 3. a) The MLE, Ai, as a function of U for various system sizes N = 50, 100, 200, 500 
(the vertical line shows the critical energy Uc = 3/4). b) MLE multiplied by the factor 
vs. U. c) MLE muhiphed by the factor A^^^ vs. U. 

Since the KS entropy is the sum of all positive Lyapunov exponents, inheriting some 
features of the LS, we expect that the scaling should be revealed more easily. This 
is indeed the shown in Fig. HI where we plot the KS entropy density Sks/N 

vs. U for several system sizes (panel a) and the same quantity multiplied either by 
N^^^ (panel b) or by N^^^ (panel c). From panel a) we observe that the KS entropy 
shows a general trend to decrease for all energies, as the system size increases. This is 
more evident than for the MLE. Panel b) shows that this size dependence is almost 
completely eliminated for U < Uc, hence the A^~^/^ factor captures some essential 
feature of the scaling for all energies. However, in the U > Uc range we still observe 
a significant size dependence. This latter is almost completely eliminated in panel 
c), where we multiply the KS entropy by N^/^, a phenomenological factor for which 
we have no theoretical justification. The low energy KS entropy scaling is however 
spoiled by this rescaling, showing that a global scaling, valid for all energies, does 
not seem to exist. It is also clear that the KS entropy peak is shifted at energies 
larger than Uc, at variance with the behavior of the MLE. 



N-50 
M=IUO 
N=200 
N-500 





Fig. 4. a) KS entropy density vs. energy U for the same data and system sizes as in Fig. 
El b) KS entropy density vs. U multiplied by the factor N^/'^ . c) KS entropy density vs. 
U multiplied by the factor iV^/^. 



The dependence of t he MLE on system size will be investigated in more detail in 



Leoncini et al.l (120101 ). going to number of particles of the order of = 10^. 
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5 The GALI indices 



In this Section we will briefly introd uce the deflnition of the Generalized Alignment 
Indices (GALI) f lSkokos et al.l . 120071 ) and we will comment on their physical meaning 
and their relation with Lyapunov exponents. We will also present the application 
of the method to the HMF model with a limited number of particles, to show its 
efficiency for the fast detection of chaotic orbits. 

The GALI index of order p (GALIp) is determined through the evolution of 2 < 
p < 2N initially linearly independent deviation vectors Wj(0),2 = 1,2,..., p. Its 
calculation is therefore strongly related to the one of the LS. The evolved deviation 
vectors Wj(t) are normalized at given time intervals in orde r to avoid ov e rflows , 
but their directions are kept unchanged. Then, according to ISkokos et al.l (l2007l ). 
GALIp is deflned to be the volume of the p-parallelogram having as edges the p unit 
deviation vectors Wj(t) = Wj(t)/||wj(t)||, i = 1,2, ...,p 



GALIp(t) =11 wi(t) A W2(t) A ... A Wp(t) II . (18) 

From this deflnition it is evident that if at least two of the deviation vectors be- 
come linearly dependent during dynamical evolution, the wedge product in Eq. fll8p 
becomes zero and the GALIp vanishes. 

In the case of a chaotic orbit, all deviation vectors tend to become linearly dependent, 
aligning along the eigenvector corresponding to the ma ximal Lyapuiiov exp onent and 



GALIp tends exponentially to zero following the law flSkokos et al.l . |200 



m 



GALIp (t) 



~ e 



-[(Ai-A2)+(Ai-A3)+...+(Ai-Ap)]f 



(19) 



where the Aj are the Lyapunov exponents (or, better, their flnite time approxima- 
tions). 

In the case of regular motion, all deviation vectors move on the A^-dimensional tan- 
gent space of a torus, on which the motion is quasiper iodic. Thus, if one starts with 
p < N deviation vectors, these will remain linearly independent on this tangent 
space, since there is no particular reason for them to become aligned. As a con- 
sequence, GALIp remains in this case practically constant for all p < N. On the 
other hand, for p > N, GALIp tends to zero, sinc e some deviation vector w il l even - 
tually becorne line arly dependent. According to IChristodoulidi &: Bountid (120061 ): 
Skokos et al.l (120081 ) . the GALI indices associated with a quasiperiodic orbit lying 
on a m-dimensional torus (with m < N) behave as follows 



GALIp (t) 



constant, if 2 < p < m 
j^^, a m < p < 2N — m 
ii2N -m<p<2N. 



(20) 
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When m = N, GALIp remains constant for 2 < p < N and decreases to zero as 
~ 1 for < p < 2N. An efficient way to calculate GALIp is by multiply- 

ing the singular values Zi,i = l,...,p, computed using the singular value decom- 
position procedure of the matr i x formed by the de viation vectors = l,...,p 



flAntonopoulos fc Bountid . l2006l : ISkokos et al.l . [20081) 



GALIp = n 



Zi. 



(21) 



i=l 



The method has been app l ied su ccessfully to several Hamiltonian s ystems like the 



FPU lattice (jSkokos et al.l . l2008l ) and to coupled symplectic maps (IBountis et al. 



20091 ). It has been shown that it efficiently detects not only regular and chaotic 



motion but also the dimensionality of the tori on which the regular trajectory lies. 

We begin by analyzing the GALI indices for systems with a few particles and we show 
examples of computations performed for the HMF model with N = 2 (integrable 
case) and A^ = 3. 

For A^ = 2 we have access to GALl2,3,4, since in this case the dimension of the 
phase space is 2A^ = 4. Besides energy, also total momentum is conserved; hence 
the system is integrable and the motion is expected to be regular for all energies. 
Following the definition of the GALI indices for regular motion, Eq. (l20l) . we should 
expect 

GALUt) oc const., GALUt) oc ^, GALl4(t) oc ^. (22) 

Indeed, choosing an initial condition with U = 0.35 and integrating it up to t = 10^, 
we see in panel a) of Fig. [5] that it produces exactly the predicted time evolution for 
the GALl2,3,4. 

For A^ = 3 we detect both regular motion (Fig. [Sb for U = 0.07 and Fig. (SJi for 
U = 1.7) and chaotic motion (Fig. |St for t/ = 0.51). In the regular case GALl2,3 
stay constant, implying a motion that lies on a 3 dimensional torus while GALI45 g 
decay following the power laws 

GALl4(t) oc ^, GALl5(t) oc ^, GALl6(t) oc ^. (23) 



Chaotic orbits of 3-dimensional Hamiltonian systems generally have two positive 
Lyapunov exponents. However, in the HMF model also total momentum is con- 
served and, therefore, only one positive exponent Ai is present in the spectrum. 
Furthermore, due to the symmetry of the Lyapunov spectrum Xq = — Ai. Therefore, 
according to formula (fT9|) . the exponential decay in time of the GALI indices is 
controlled only by Ai, as shown in Fig. [St. 
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"t Log(t) 



Fig. 5. GALI indices vs. time for regular and chaotic orbits, a) HMF model with N = 2 
and energy U = 0.35: integrable system for which GALl3^4 decay to zero with the power 
law predicted theoretically, b) HMF model with = 3 and energy U = 0.07 for which 
GALl2,3 are constant: the motion lies on a 3 dimensional torus, c) HMF model with N = 3 
and energy U = 0.51: chaotic motion is detected, i.e. all GALI's decay exponentially to 
zero. The slopes are given by the Lyapunov exponents, d) HMF model with N = 3 and 
energy U = 1.7: the motion is regular and lies again on a 3-dimensional torus. 

The most important advantage of the GALI method rehes on the fact that, for 
chaotic orbits, the GALI indices decay to zero exponentially fast. Practically, an 
orbit can be labeled as chaotic when some GALIp (where p is conveniently chosen) 
is found to be smaller than a given threshold. This often happens before the MLE 
stabilizes on an average positive value, making the GALIp test more efficient in 
determining chaotic properties of an orbit than the measurement of the MLE. We 
will give an explicit example of this property of the GALI indices in the next Section. 



6 Application of the GALI method to the HMF model with large A^: 
the transition from weak to strong chaos 



In this Section we will present the application of the GALI method to the HMF 
model, as a tool to reveal the fraction of chaotic orbits on a constant energy surface 
with many degrees of freedom. The dependence of such a fraction on the number 
of degrees of freedoi n has been re c ently analyzed for fully or partially connected 
symplectic maps in iLaveder et al.l ( l2008l ). with an application also to the HMF 
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model. 



As we have discussed in the previous Section, the GALI method can be faster in 
reveahng the chaotic or regular nature of a given orbit than the calculation of the 
MLE. In order to perform a systematic analysis, it is necessary to fix some criteria. 
The first important choice is the value of p for the GALIp. In systems with many 
degrees of freedom we have access to high values of p, but the calculation of the 
corresponding GALIp would be extremely heavy, since we would have to perform the 
singular value decomposition for large matrices. Therefore, we are limited ourselves 
to small values of p. We have avoided to use GALI2, since it has been shown in 



Skokos et al.l (120071 ) that, for some chaotic orbits of many degrees of freedom systems, 
the two largest Lyapunov exponents can take very close values. In such chaotic cases, 
the GALI2, instead of decreasing exponentially, stays constant, detecting a "false" 
regular orbit. For the HMF model, in the energy range < U < 0.4, we have 
checked the different time evolution of GALI345 and we have finally decided to 
concentrate our attention on the behavior of GALI3. Fixing the total integration 
time to t = 3 X 10^, we have labeled an orbit as regular if GALI3 stays above 
10~^ over the whole time lapse. We instead label an orbit as chaotic if GALI3 goes 
below 10~^^ within the total integration time. Orbits not satisfying any of these two 
criteria are analyzed on a longer time span to assess their chaotic or regular nature; 
however, in the simulations we have performed very few such cases happened. We 
have decided to fix the number of initial conditions to 1000 and we have performed 
simulations with increasing number of particles: = 100, 1000, 5000. 

In Fig. [6] we show the fraction of chaotic orbits in % as a function of U. A sharp 
increase of this fraction is observed around the energy Ut ~ 0.2. A more quantitative 
fit to the data for = 1000 with the function /(x) = a * [1 + tanh(7 * {U — Ut))] 
gives the more precise estimate Ut = 0.217... (the other fitting parameters being 
a = 50 and 7 = 30). A significant shift of the transition towards smaller values 
of U is observed as A^ is increased. However, it is likely that the curve will reach 
an asymptotic limit as A^ is increased (currently not accessible for our computers), 
since the shift observed when passing from A^ = 100 to A^ = 1000 (a factor of 10) 
is of the same order as the one found when further multiplying A^ by a factor of 5. 
Moreover, the fraction of chaotic orbits is below 1% for f/ < 0.15 and around 99% 
for U > 0.25 for all simulated system sizes, showing that a sharp transition from 
a weakly chaotic to a strongly chaotic phase-space takes place in the energy range 
[0.15,0.25]. 



It has been suggested by lAntoni fc Ruffd f ll995l ) that, around energy 0.2, the phase 



(j) of the order parameter M of the HMF model begins to be strongly time depen- 
dent, while for smaller energies it is essentially time independent. This determines 
the untrapping of some particle trajectories from the nonlinear resonance in the pen- 
dulum phase-space associated to the model. Since we here find a significant growth 
of the fraction of chaotic orbits just in this energy range, we wanted to verify the 
correspondence between the two phenomena. 
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u 

Fig. 6. Fraction of chaotic orbits vs. energy U for a different number of particles 
N = 100,1000,5000. We also show a fit to the points with N = 1000 particles (full 
line). The transition energy from weak to strong chaos is around Ut ~ 0.2. 



In Fig. [7] (upper row) we plot the time evolution of GALI3 for three different energies: 
U = 0.075, 0.2, 0.4. At low energy (panel a) the GALI3 remains constant, revealing a 
regular orbit, while at high energy (panel c) it decreases exponentially fast, detecting 
a chaotic orbit. In the energy region around Ut, part of the orbits are chaotic and 
part are regular. We here show (panel b), two orbits that, according to our criteria, 
are labeled as chaotic (thick full line) and regular (thin full line). In the inset of 
panel b) we plot the time evolution of the MLE for these two orbits: the MLE of 
the chaotic orbit is larger that the one of the "regular" orbit, but on this time scale, 
they are still far from convergence. This proves that the GALI3 is able to detect the 
regular and chaotic nature of the orbit faster than the MLE. 

Phase motion is absent for the regular orbit at low energy, see Fig. [TJi, while it's 
clearly present for the chaotic orbit at the highest energy in Fig. [Tf. At the energy 
U = 0.2 (Fig. [7^) both the chaotic and the regular orbit display phase motion but, 
surprisingly, the regular orbit shows a "ballistic" motion of the phase. It's unclear 
why this happens; a possible explanation is that, for this orbit, untrapped particles 
remain out of resonance for a longer time. This would imply less chaoticity because 
some particles remain far from the most chaotic part of the phase-space, which is 
around the chaotic layer of the separatrix. 

Although the understanding of the transition from weak to strong chaos appearing 
around the energy Ut is far from being satisfactory, its numerical evidence is clear. 
Moreover, this threshold energy separates an energy region, < U < Ut, where the 
A^~^/^ scaling of the MLE and of the LS is clearly found from a higher energy region, 
Ut < U < Uc, where this scaling is not found with the same evidence. Indeed, in this 
intermediate energy range, the MLE hardly converges to zero as the system size is 
increased, as shown in Fig. [3^. 
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Fig. 7. GALI3 vs. time (upper row) and phase of the order parameter M (see Eq. [2]) vs. 
time (lower row) at increasing energies for N = 1000. a) GALI3 remains constant in time 
for a regular orbit at U = 0.075. b) GALI3 for a chaotic orbit (thick full line) and for a 
less chaotic one, classified as regular according to our criterion (thin full line), at U = 0.2. 
In the inset we show the corresponding time evolution of the MLE. c) GALI3 decays 
exponentially for a typical chaotic orbit at U = 0.4. d) Phase (p of the order parameter vs. 
time for U = 0.075 computed for the orbit of panel a), e) Phase (j) of the order parameter 
vs. time for U = 0.2 computed for the two orbits of panel b). f) Phase (j) of the order 
parameter vs. time for U = 0.4 computed for the chaotic orbit of panel c). 

7 Conclusions 



Although the chaotic properties of the Hamiltonian Mean Field (HMF) model have 
been studied in several papers, a general understanding of the behavior of the Lya- 
punov exponents as the system size is increased and for different energies is still 
lacking. 

We have here shown that, rather than concentrating on the Maximal Lyapunov Ex- 
ponent (MLE), it is sometimes preferable to determine the full Lyapunov Spectrum 
(LS). This is certainly more challenging from the numerical point of view, but can 
be rewarding and reserve some pleasant surprises. We here show, for instance, that 
the convergence of the LS to its asymptotic, large N, shape is more rapid than the 
one of the MLE. 



Moreover, b y the study of the LS, the general scaling law with N~^^^, originally 
proposed by iLatora et all (119981 ). is here found for moderate system sizes. Although 
this law was originally proposed for the high energy phase of the HMF, above the 
critical energy t/c = 3/4, we show in this paper that it is also valid in a wide range 



15 



of low energies, below Ut ~ 0.2. 



With the aim of investigating the physical meanin g of this threshold e nergy Ut, we 
have used the newly proposed GALI method by ISkokos et al.l (120071 ) in order to 
determine the fraction of chaotic orbits in the phase-space of the HMF model. The 
method reveals that this fraction sharply changes from few percents to almost 100% 
in an energy range around Ut, showing a transition from weak to strong chaos in 
this energy region. It is for energies < U < Ut that the N~^^^ law is best verified. 
On the contrary, when Ut < U < Uc, the scaling law is numerically unclear and the 
MLE shows a slow decrease to zero as the system size increases, which could even 
lead to suspect that the MLE remains positive in the N ^ oo limit. 

Above Uc one observes a general trend to decrease to zero of the full LS. The scaling 
is precocious if one studies full the LS and is again well fitted by the N~^^^ law. 
The MLE also decreases to zero, but with an intermediate scaling with N~^^^. It is 
only for much larger system sizes, of the order of = 10^, that one finally finds the 
N-'^/^ law also for the MLE. 



On the theor etical side, the only existing approach valid in principle for all energies 
( iFirpoi Il998l ) predicts a strictly positive MLE for < f/ < Uc, while the MLE 
should vanish as A^~^/^ for U > Uc- This is not what we find for < U < Ut, 
demanding for an improvement of this theoret ical approach. In the high energy 



regio n U > Uc, a random ma t rix approximation (ILatora et al.l. Il998l: iFirpo fc Ruffo 



2001 



Anteneodo fc Vallejosl . 1200 ll : IVallejos Sz Anteneodd . l2002l : lAnteneodo et al. 



2003) predicts t h e A^ scaling of the MLE , which has been checked nume rically 



bv iLatora et al.l f ll998l ): lAnteneodo fc TsallisI f ll998l ): lAnteneodo et al.l ( 120031 ). It is 
hard to believe that such an approach could be used for the energy range < U < Ut, 
where the phase-space is mostly occupied by regular orbits, as the GALI method 
allowed us to assess. The understanding of why the A^~^/^ scaling is valid in this 
energy region remains a challenge for future research. 
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